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Abstract. We explore the regular or chaotic nature of orbits of stars mov¬ 
ing in the meridional {R, z) plane of an axially symmetric time-dependent disk 
galaxy model with a central, spherically symmetric nucleus. In particular, mass 
is linearly transported from the disk to the galactic nucleus, in order to mimic, 
in a way, the case of self-consistent interactions of an actual N-body simulation. 
We thus try to unveil the influence of this mass transportation on the different 
families of orbits of stars by monitoring how the percentage of chaotic orbits, as 
well as the percentages of orbits of the main regular resonant families, evolve as 
the galaxy develops a dense and massive nucleus in its core. The SALI method 
is applied to samples of orbits in order to distinguish safely between ordered 
and chaotic motion. In addition, a method based on the concept of spectral 
dynamics is used for identifying the various families of regular orbits and also 
for recognizing the secondary resonances that bifurcate from them. Our compu¬ 
tations strongly suggest that the amount of the observed chaos is substantially 
increased as the nucleus becomes more massive. Eurthermore, extensive numer¬ 
ical calculations indicate that there are orbits which change their nature from 
regular to chaotic and vice versa and also orbits which maintain their orbital 
character during the galactic evolution. The present outcomes are compared to 
earlier related work. 

Key words: galaxies: kinematics and dynamics - structure - chaos 
1. INTRODUCTION 


A convenient way for studying autonomous Hamiltonian systems is by choos¬ 
ing some values of the total orbital energy, for which the position as well as the 
extent of ordered and chaotic domains are time-independent (TI) and can be accu¬ 
rately determined by a variety of dynamical methods, especially in the case of low 
degrees of freedom. The dynamical structure however, in such time-independent 
Hamiltonian systems is usually very complex due to the presence of “weak” and 
“strong” chaos. In addition, the motion of test particles takes place on invariant 
tori of multiple dimensions and it can also display surprising localization proper¬ 
ties in both configuration and phase space. It is therefore natural to assume that 
time-dependent (TD) Hamiltonian systems are expected to exhibit a much more 
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complicated nature since all the above-mentioned attributes evolve with time due 
to the absence of any kind of integrals of motion (i.e., the total energy integral) 
which are valid only in TI systems. Indeed, it is true that in TI systems orbits 
maintain their orbital character: if for example, an orbit is initially regular it will 
always remain so, while if it is chaotic it might get trapped near the boundaries 
of stability regions for relatively long time intervals (this phenomenon is known 
as “sticky motion”) but it will certainty never relinquish its chaoticity. On the 
other hand, in TD Hamiltonians sudden transitions form regularity to chaotic¬ 
ity and vice versa is a common behavior for individual trajectories during their 
time-evolution. 

The determination of ordered and chaotic properties of motion in time-dependent 
galactic as well as cosmological models constitutes an extended research area in 
the field on non-linear dynamics. In a previous work, Caranicolas & Papadopoulos 
(2003) used a simple analytic time-dependent model in order to study the tran¬ 
sition from order to chaos in a galaxy, when mass is exponentially transported 
from the disk to the galactic core thus forming a dense and massive spherical 
nucleus. They found that during the galactic evolution a large portion of low an¬ 
gular momentum stars change their orbital nature from regular to chaotic. The 
investigation was continued in more detail in Zotos (2012) where it was proved 
by conducting a systematic and thorough exploration of the phase space that 
during the mass transportation stars do not change their character only from reg¬ 
ular to chaotic. In fact, there is a considerable amount of stars which maintain 
their orbital nature, while there is also a small fraction of star orbits that change 
their character from chaotic to regular. Furthermore, Papadopoulos & Caranico¬ 
las (2006) revealed the orbital behavior in a time-dependent double-barred galaxy 
model when mass is transported from the primary bar to the inner disk reporting 
that the galactic evolution significantly affects the nature of orbits. The interplay 
between ordered and chaotic motion in a time-dependent Hamiltonian describing 
a barred galaxy was examined in Manos et al. (2013), where it was observed that 
the percentage of chaos increases when a linear mass transportation from the disk 
to the bar takes place. In the same vein, Manos & Machado (2014) constructed 
an analytical time-dependent potential, modeling the gravitational potentials of 
disk, a bar and a dark matter halo, whose time-dependent parameters are derived 
from a simulation. 

Knowing whether the orbits are regular or chaotic is only the first step to¬ 
ward the understanding of the overall behavior of the disk galaxy. The second 
and beyond any doubt the most interesting step is the distribution of regular 
orbits into different families. In our study therefore, once the orbits have been 
characterized as regular or chaotic, we then further classified the regular orbits 
into different families by using a frequency analysis method (Carpintero & Aguilar 
1998; Muzzio et al. 2005). Initially, Binney & Spergel (1982, 1984) proposed a 
technique, dubbed spectral dynamics, for this particular purpose. Later on, this 
method has been extended and improved by Sidlichovsky & Nesvorny (1996) and 
Carpintero & Aguilar (1998). In the recent work of Zotos & Carpintero (2013) the 
algorithm was refined even further so it can be used to classify orbits in the merid¬ 
ional plane. In general terms, this method computes the Fourier transform of the 
coordinates of an orbit, identifies its peaks, extracts the corresponding frequencies, 
and searches for the fundamental frequencies and their possible resonances. Thus, 
we can easily identify the various families of regular orbits and also recognize the 
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secondary resonances that bifurcate from them. This technique has been success¬ 
fully applied in several previous papers (e.g., Zotos & Carpintero 2013; Caranicolas 
& Zotos 2013; Zotos & Caranicolas 2013, 2014; Zotos 2014a,b) in the field of orbit 
classification (not only regular versus chaotic, but also separating regular orbits 
into different regular families) in different galactic gravitational potentials. 

The motivation of the present work is to apply the time-dependent version of 
the disk galaxy model used in Zotos & Carpintero (2013) in order to examine the 
dynamical properties and the time-evolution of orbits as mass is transferred from 
the disk to the central nucleus. Moreover, this time-dependent Hamiltonian can, 
in a way, mimic certain realistic aspects that arise in N-body simulations. Our 
paper is organized as follows: in Section 2, we explain in detail the properties of the 
time-dependent disk galaxy model, while Section 3 is devoted to the description 
of the computational methods we used in order to determine the character as 
well as the classification of the orbits. A thorough numerical analysis regarding 
the influence of the mass transport to the percentages of the different families of 
orbits is performed in Section 4. Our paper ends with Section 5, where the main 
conclusions of our numerical investigation are presented. 

2. DESCRIPTION OF THE GALACTIC MODEL 


Our aim is to determine the interplay between ordered and chaotic motion of 
stars moving in the meridional plane of a time-dependent axially symmetric disk 
galaxy with a central spherically symmetric nucleus when the parameters of the 
model potential evolve in time. For this purpose, we use the usual cylindrical 
coordinates (R, 0, z), where z is the axis of symmetry. 

The total gravitational potential 4>(R, z) consists of two components: the nu¬ 
cleus potential 4>n and the flat disk potential 4>d. For the description of the 
spherically symmetric central nucleus, we use a Plummer potential (e.g., Binney 
& Tremaine 2008) 


^n{R,z) 


-GMn 

yi?2 + + c2 


( 1 ) 


Here G is the gravitational constant, while and Cn are the mass and the scale 
length of the nucleus, respectively. This potential has been successfully used in the 
past to model and, therefore, interpret the effects of the central mass component 
in a galaxy (see e.g., Hasan & Norman 1990; Hasan et al. 1993; Zotos 2012; Zotos 
& Carpintero 2013; Zotos 2014a). At this point, it must be emphasized that we do 
not include any relativistic effects, because the nucleus represents a bulge rather 
than a black hole or any other compact object. 

The galactic disk, on the other hand, is represented by the well-known Miyamoto- 
Nagai potential Miyamoto & Nagai (1975) 


^d{R,z) 


_ -GMj _ 

\jB? + {a-\- + z^) 


( 2 ) 


where, is the mass of the disk, a is the scale length of the disk and h corresponds 
to the disk’s scale height. 

We use a system of galactic units where the unit of length is 1 kpc, the unit of 
velocity is 10 kms“^ and G = 1. Thus, the unit of mass is 2.325 x IO^Mq, that 
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of time is 0.9778 x 10^ yr, the unit of angular momentum (per unit mass) is 10 
km kpc s“^, and the unit of energy (per unit mass) is 100 km^s“^. In these units, 
the values of the involved parameters are: Cn = 0.25, a = 3, h = 0.175, while 
Mn and Md are treated as variable parameters with time. The above-mentioned 
set of values of the parameters which are kept constant throughout the numerical 
calculations secures positive mass density everywhere and free of singularities. 

For simplicity, we assume that mass is transported from the disk to the nucleus, 
in such as way, that we have a linear increase in the mass of the nuclear core, while 
a simultaneous linear decrease in the mass of the disk takes place. The linear mass 
transportation follows the equations 


Tfn — Tfno + ^ ^5 

Md = Mdo - k t, (3) 


where Mno and Mdo are the initial values of the mass of the nucleus and the 
disk, respectively, while k > 0 is the proportionality constant which determines 
the timescale of the galactic evolution. We choose Mno = 0, Mdo = 7500 and 
k = 0.05, so at the beginning of the galactic evolution {t = 0) there is no nucleus 
formed and all the mass is concentrated in the disk. We also assume that the 
linear rate described by Eqs. (3) is slow compared to the orbital period of the 
stars and is therefore adiabatic. This is true because the mass transportation 
lasts for 10^ time units, while the orbital period is about two orders of magnitude 
smaller. Furthermore, the mass transportation and the galactic evolution stops 
after 10^ time units when the final value of the mass of the nucleus and that of 
the disk is Mn = 500 and Md = 7000, respectively. The particular final values of 
the parameters were chosen with a Milky Way-type galaxy in mind (e.g., Allen 
& Santillan 1991). It is well known, that such mass transportation mechanisms 
are usually met in the central regions of active galaxies and they are responsible 
for the enormous luminosity of the quasars which are hosted in the active galactic 
nuclei of the galaxies (see e.g., Collin & Zahn 1999). This is the main reason why 
we refer to the central mass concentration as a “nucleus” rather than as bulge (see 
also Zotos 2012). 

The circular velocity in the galactic plane where z = 0 (see e.g., Zotos 2011) 
is undoubtedly one of the most important physical quantities in disk galaxies and 
can be calculated as 


0{R) = 



d^{R,z) 

dR 


z=0 


( 4 ) 


In Fig. 1 we present a plot of 0{R) at the beginning (t = 0) and at the end 
(t = 10000) of the galactic evolution. The red line corresponds to the case where 
Mn = 0 and Md = 7500, that is when the central nucleus is absent, while the green 
line corresponds to the case where Mn = 500 and Md = 7000 and a fully developed 
dense, massive nucleus lies in the galactic core. It is seen that at relatively low 
galactocentric distances {R < b kpc) the pattern of the two curves is completely 
different, while at large distances from the galactic center {R > b kpc), on the 
other hand, both curves fully coincide. We also observe the characteristic local 
minimum {at R 1 kpc) of the rotation curve when the massive nucleus is present, 
which appears when fitting observed data to a galactic model (e.g., Gomez et al. 
2010; Irrgang et al. 2013). 
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Fig. 1. A plot of the total circular velocity of our galactic model when Mn = 0 and 
Md = 7500 (red) and when Mn = 500 and Md = 7000 (green). 


Exploiting the fact that the L; 2 -component of the total angular momentum is 
conserved because the gravitational potential ^{R,z) is axially symmetric, orbits 
can be described by means of the effective potential 

^,ii{R,z) = ^R,z) + fT, ( 5 ) 


In this case, the basic equations of motion on the meridional plane take the 
form 


dR ’ ^ dz 


(6) 


while the equations governing the evolution of a deviation vector w = {6R, 6z^ 6R, Sz), 
which joins the corresponding phase space points of two initially nearby orbits, 
needed for the calculation of the standard indicators of chaos (the SALI in our 
case), are given by the following variational equations 


(SR) 

(SR) 

(Sz) 


SR, (Sz) = Sz, 

dR^ dRdz 

dzdR dz‘^ 


Sz, 

Sz. 


( 7 ) 


Consequently, the corresponding Hamiltonian to the effective potential given 
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in Eq. (5) can be written as 

H=\{R^+z‘^)+^es{R,z) = E, (8) 

where R and i are momenta per unit mass, conjugate to R and z, respectively, 
while E is the numerical value of the Hamiltonian, which is conserved only in 
the time-independent model. Therefore, an orbit is restricted to the area in the 
meridional plane satisfying E > while all other regions are forbidden to the 
star. 

3. COMPUTATIONAL METHODS 

In order to examine the orbital dynamics (regularity or chaoticity) of the galaxy 
model, we need to establish some samples of initial conditions of orbits. The best 
approach, undoubtedly, would be to extract these samples of orbits from the dis¬ 
tribution function of the model potential. Unfortunately, this is not available, so 
we followed another course of action. To determine the character of the orbits in 
our model, we chose dense grids of initial conditions in the phase (i?, R) space reg¬ 
ularly distributed in the area allowed by the value of the value of the Hamiltonian 
E. At this point we have to point out that the value of the energy controls the 
grid size and particularly the i^max which is the maximum possible value of the 
R coordinate. On this basis, we chose as initial energy level the value Eq = —480 
which yields i^max — 15 kpc. Moreover, the angular momentum of all orbits is 
Lz = 10 and remains constant throughout. Since the Hamiltonian of our model 
is time-dependent, its value is not conserved during the galactic evolution. To 
determine the change in the value of the Hamiltonian, i.e. dE = {E — Eq)/Eq, we 
chose 10000 random initial conditions in the phase space with Eq = —480 and we 
recorded the value of the Hamiltonian of every orbit during the galactic evolution 
until t = 10^ time units when the mass transportation stops. Our results are 
presented in Eig. 2 where we observe the evolution of the value of Hamiltonian as 
a function of time. It was found that for all integrated orbits their energy is lin¬ 
early reduced with time within the gray-shaded area of the plot. Eurthermore, the 
change of the value of the Hamiltonian ranges from about 3.5% to 4.8% and the 
black, dashed lines in Eig. 2 indicate the corresponding minimum and maximum 
linear trends. 

Taking into account that the range of values regarding the change of the Hamil¬ 
tonian is relatively small, we decided to adopt an average linear change for all orbits 
in our model. Our calculations show that the average Hamiltonian change, shown 
in blue in Eig. 2, follows the linear rule 


E = Eo-qt, (9) 

where Eq = —480, while q = 0.002 so that E = —500 at t = 10^ time units. 
Relation (9) is a powerful tool because it allows us to know the exact value of 
the Hamiltonian at every time step of the galactic evolution. Now we can use 
this valuable information in order to explore the orbital structure of the galaxy at 
various time points during the mass transportation. In particular, since we have 
equations describing the time-evolution of Mn, Md and E, we can pick some time 
points ti, freeze the mass transportation and use the time-independent model in 
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Fig. 2. Evolution of the value of the time-dependent Hamiltonian as a function of 
time for 10000 orbits initiated at the phase {R,R) space with Eq = —480. The blue 
line indicates the average Hamiltonian change, while the black, dashed lines delimit the 
boundaries of the energy variation dE. 

order to integrate and classify the orbits, following the numerical approach used 
in Manos et ah (2013). Therefore, for each set of values of M^{ti) and Md(ti), we 
define a dense grid of initial conditions in the phase space regularly distributed in 
the area allowed by the corresponding value of the Hamiltonian E(ti). The step 
separation of the initial conditions along the R and R axes (or, in other words, 
the density of the grid) was controlled in such a way that always there are at 
least 50000 orbits to be integrated and classified. The grids of initial conditions 
of orbits whose properties will be examined are defined as follows: we consider 
orbits with initial conditions (Rq^Rq) with zq = 0, while the initial value of zq 
is obtained from the Hamiltonian (8). For each initial condition, we numerically 
integrated the equations of motion (6) as well as the variational equations (7) with 
a double precision Bulirsch-Stoer FORTRAN 77 algorithm (e.g.. Press et al. 1992) 
with a small time step of the order of 10“^, which is sufficient enough for the 
desired accuracy of our computations (i.e., our results practically do not change 
by halving the time step). In all cases, the energy integral^ (Eq. 8) was conserved 
better than one part in 10“^^, although for most orbits it was better than one part 
in 10“^^. 

For determining the regular or chaotic nature of orbits, we chose the SALI 


^Remember that these orbits are integrated in the frozen time-independent model, where the 
energy integral is valid. 
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indicator Skokos (2001). The time-evolution of SALI strongly depends on the 
nature of the computed orbit since when an orbit is regular the SALI exhibits 
small fluctuations around non-zero values, while on the other hand, in the case 
of a chaotic orbit, the SALI after a small transient period tends exponentially to 
zero approaching the limit of the accuracy of the computer (10“^^). Therefore, 
the particular time-evolution of the SALI allow us to distinguish fast and safely 
between regular and chaotic motion. Nevertheless, we have to define a specific 
numerical threshold value for determining the transition from regularity to chaos. 
After conducting extensive numerical experiments, integrating many sets of orbits, 
we conclude that a safe threshold value for the SALI, taking into account the total 
integration time of 10^ time units, is the value 10“^. In order to decide whether 
an orbit is regular or chaotic, one may follow the usual method according to which 
we check, after a certain and predefined time interval of numerical integration, if 
the value of SALI has become less than the established threshold value. Therefore, 
if SALI < 10 ^ then the orbit is chaotic, while if SALI > 10 ^ then the orbit is 
regular. However, depending on the particular location of each orbit, this threshold 
value can be reached more or less quickly, as there are phenomena that can hold 
off the final classification of an orbit (i.e., there are special orbits called “sticky” 
orbits, which behave regularly for long time periods before they finally drift away 
from the boundaries of regular regions and start to wander in the chaotic domain, 
thus revealing their true chaotic nature fully. For the computation of SALI we 
used the LP-VI code Carpintero et al. (2014), a fully operational routine which 
efficiently computes a suite of many chaos indicators for dynamical systems in any 
number of dimensions. 

Each orbit in the frozen time-independent Hamiltonian was numerically inte¬ 
grated for a time interval of 10^ time units (10^^ yr), which corresponds to a time 
span of the order of hundreds of orbital periods. The particular choice of the total 
integration time is an element of great importance, especially in the case of the 
sticky orbits. A sticky orbit could be easily misclassified as regular by any chaos 
indicator^, if the total integration interval is too small, so that the orbit does not 
have enough time to reveal its true chaotic character. Thus, all the initial con¬ 
ditions of the orbits of a given grid were integrated, as we already said, for 10^ 
time units, thus avoiding sticky orbits with a stickiness at least of the order of 
10^ Hubble times. All the sticky orbits that do not show any signs of chaoticity 
for 10^ time units are counted as regular orbits since such vast sticky periods are 
completely out of the scope of our research. 

A vital clarification regarding the nomenclature of the orbits should be made 
before closing this Section. All orbits of an axially symmetric potential are in 
fact three-dimensional (3D) loop orbits, i.e., orbits that always rotate around the 
axis of symmetry in the same direction. However, in dealing with the meridional 
plane, the rotational motion is lost, so the path that the orbit follows onto this 
plane can take any shape, depending on the nature of the orbit. Following the same 
approach of the previous papers of this series, we characterize an orbit according 
to its behavior in the meridional plane. If, for example, an orbit is a rosette lying 
in the equatorial plane of the axisymmetric potential, it will be a linear orbit in 


^ Generally, dynamical methods are broadly split into two types: (i) those based on the evo¬ 
lution of sets of deviation vectors to characterize an orbit and (ii) those based on the frequencies 
of the orbits that extract information about the nature of motion only through the basic orbital 
elements without the use of deviation vectors. 
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Table 1. Types and initial conditions of the orbits shown in Fig. 3(a-h). In all cases, 
zq = 0, zq is found from the Hamiltonian, Eq. (8), while Tper is the period of the resonant 
parent periodic orbits. 


Figure 

Type 

Ro 

Ro 

Tper 

3a 

box 

1.20000000 

0.00000000 

- 

3b 

2:1 banana 

6.94608753 

0.00000000 

1.82936813 

3c 

1:1 linear 

3.48802471 

38.08701776 

1.52281338 

3d 

2:3 boxlet 

14.12356085 

0.00000000 

3.07990781 

3e 

4:3 boxlet 

12.62569659 

0.00000000 

5.95238344 

3f 

6:5 boxlet 

13.27429077 

0.00000000 

9.51868683 

3g 

4:5 boxlet 

14.30689514 

0.00000000 

6.47237905 

3h 

chaotic 

0.12000000 

0.00000000 

- 


the meridional plane, a tube orbit it will be a 2:1 orbit, etc. We should emphasize 
that we use the term “box orbit” for an orbit that conserves circulation, but this 
refers only to the circulation provided by the meridional plane itself. Because of 
the their boxlike shape in the meridional plane, such orbits were originally called 
“boxes” (e.g., Ollongren 1962), even though their three-dimensional shapes are 
more similar to doughnuts (more details can be found in the review of Merritt 
1999). Nevertheless, we kept this formalism to maintain continuity with all the 
previous papers of this series. 

4. NUMERICAL RESULTS 


This Section contains the main numerical results of both the time-independent 
and the time-evolving model. In particular, we numerically integrate several sets 
of orbits, in an attempt to determine the regular or chaotic nature of motion 
of stars. First we begin with the time-independent model and use the sets of 
initial conditions described in the previous section to construct the respective 
grids, always adopting values inside the limiting curve. In all cases, the initial 
value of the energy was set equal to —480, while the angular momentum of the 
orbits is Lz = 10. To study how the mass transportation from the disk to the 
central nucleus influences the level of chaos, we choose representative time points 
of the galactic evolution such as U = {0, 2000,4000,..., 10000}. Once the exact 
values of the parameters Mn, Md and E are known through Eqs. (3) and (9), we 
compute a set of initial conditions as described in Section 3 and we integrate the 
corresponding orbits computing the SALI of the orbits and then classifying regular 
orbits into different families. 

Our numerical investigation suggests that in our galaxy model there are eight 
basic types of orbits: (i) chaotic orbits; (ii) box orbits; (hi) 1:1 linear orbits; (iv) 
2:1 banana-type orbits; (v) 2:3 fish-type orbits; (vi) 4:3 resonant orbits; (vii) 6:5 
resonant orbit and (viii) orbits with other secondary resonances (i.e., all resonant 
orbits not included in the former categories). It turns out that for each resonant 
family included in the “other category” the corresponding percentage is less than 
1% in all cases, and therefore their contribution to the overall orbital structure of 
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the galaxy is practically insignificant. The n : m notation^ we use for the regular 
orbits is according to Carpintero & Aguilar (1998) and Zotos & Carpintero (2013), 
where the ratio of those integers corresponds to the ratio of the main frequencies 
of the orbit, where main frequency is the frequency of greatest amplitude in each 
coordinate. Main amplitudes, when having a rational ratio, define the resonances 
of an orbit. In Fig. 3(a-h) we present examples of each of the basic types of 
regular orbits, plus an example of a chaotic one. In all cases, we set Mn = 500 
and E = —500 (except for the 4:5 and 6:5 resonant orbits, where Mn = 100 and 
E = —484). The orbits shown in Figs. 3a and 3h were computed until t = 100 
time units, while all the parent periodic orbits Zotos (2013) were computed until 
one period has completed. The outermost black curve circumscribing each orbit is 
the limiting Zero Velocity Curve (ZVC) in the meridional plane which is defined 
as 4>eff(i?, 2 :) = E. Table 1 shows the types and the initial conditions for each of 
the depicted orbits; for the resonant cases, the initial conditions and the period 
Tper correspond to the parent^ periodic orbits. 

We would like to point out at this point, that the 1:1 resonance is usually the 
hallmark of the loop orbits and both coordinates oscillate with the same frequency 
in their main motion. Their mother orbit is a closed loop orbit. Moreover, when 
the oscillations are in phase, the 1:1 orbit degenerates into a linear orbit (the 
same as in Lissajous figures made with two oscillators). In our meridional plane, 
however, 1:1 orbits do not have the shape of a loop. In fact, their mother orbit 
is linear (e.g.. Fig. 3c), and thus they do not have a hollow (in the meridional 
plane), but fill a region around the linear mother, always oscillating along the R 
and z directions with the same frequency. We designate these orbits “1:1 linear 
open orbits” to differentiate them from true meridional plane loop orbits, which 
have a hollow and also always rotate in the same direction. 

In the following Figs. 4(a-f) we present six grids of initial conditions (i?o, Ro) of 
orbits that we have classified for different values of the final mass of the nucleus Mn 
in the frozen time-independent model. These color-coded grids of initial conditions 
are equivalent to the classical Poincare Surfaces of Section (PSS) and allow us to 
determine what types of orbits occupy specific areas in the phase (i?, R) plane and 
also to follow the changes that orbits undergo. The outermost black thick curve 
is the limiting curve which is defined as 

^R^ + ^^s{R,z = 0) = E. ( 10 ) 

In Fig. 4a which corresponds to time point t = 0 where the central nucleus is absent 
(Mn = 0) we observe that the entire phase plane is covered by initial conditions 
of regular orbits, while chaotic motion, if any, is negligible. In particular, initial 
conditions of box orbits occupy the vast majority of the phase space, a well-formed 
island of 2:1 resonant orbits is present at the center of the grid, while there are 
also several smaller stability islands of other resonant orbits embedded in the 
extended box area. The structure of the phase plane however, changes drastically 

^ A n : m resonant orbit would be represented by m distinct islands of invariant curves in the 
(R, R) phase plane and n distinct islands of invariant curves in the (z, z) surface of section. 

^ For every orbital family there is a parent (or mother) periodic orbit, that is, an orbit that 
describes a closed figure. Perturbing the initial conditions which define the exact position of 
a periodic orbit we generate quasi-periodic orbits that belong to the same orbital family and 
librate around their closed parent periodic orbit. 
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Fig. 4. Orbital structure of the phase {R, R) plane of our galaxy model at several 
time points of the galactic evolution. 

in Fig. 4b at t = 2000 time units when = 100. It is seen that the box 
area is significantly reduced and especially at the outer parts of the phase plane 
there is a strong presence of stability islands of resonant orbits surrounded by 
a thin chaotic layer. It is interesting to note that inside this chaotic layer one 
may identify several sets of tiny stability islands corresponding to higher resonant 
orbits which in our classification belong to the “other” category. The area on the 
phase plane occupied by these secondary higher resonant orbits is considerably 
confined when t = 4000 time units and Mn = 200. Indeed we see in Fig. 4c that 
the corresponding stability islands are hardly visible, while at the same time the 
chaotic zone is amplified. The same pattern continues in Fig. 4d for t = 6000 
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time units and Mn = 300 where the chaotic area is more prominent, while some 
thin filaments of initial conditions of secondary resonant inside the box region are 
observed. The mass transport is still in progress and for t = 8000 time units which 
corresponds to = 400 there is a major difference in the phase space, that is 
the absence of the 2:3 stability islands. Additional numerical calculations reveal 
that for Mn = 400 the 2:3 resonance is unstable. This means that the periodic 
point of the 2:3 resonance is indeed present in Fig. 4e, although evidently deeply 
buried in the chaotic domain. Finally, in Fig. 4f, where t = 10000 time units and 
Mn = 500, we observe that the 2:3 stability islands reappear in the phase space. 

Looking carefully the color-coded grids presented in Figs.4(a-f) we can distin¬ 
guish the location of the seven main types of regular orbits discussed earlier: (i) 
2:1 banana-type orbits located to the central region of the grid; (ii) box orbits 
situated mainly outside of the 2:1 resonant orbits; (hi) 1:1 open linear orbits form 
the elongated island of initial conditions; (iv) 2:3 fish-type resonant orbits form 
the set of three small islands^ at the outer parts of the grid; (v) 4:3 resonant orbits 
form the chain of three islands; (vi) 6:5 resonant orbits producing five small stabil¬ 
ity islands inside the chaotic region and (vii) other types of resonances producing 
extremely small islands embedded both in the chaotic and box areas. Chaos on 
the other hand, is mainly confined only at the outer parts of the phase plane. It 
is evident, that as the central nucleus becomes more and more massive gaining 
mass from the galactic disk, the amount of chaos in the phase space increases. It 
should also be mentioned that the allowed radial velocity R of the stars near the 
center of the galaxy is increasing during the galactic evolution, where the mass of 
the nucleus also increases. 

The time-evolution of the percentages of the chaotic and all types of regular 
orbits as a function of the final mass of the central nucleus Mn is presented in Fig. 
5. One may observe that at the beginning of the galactic evolution (t = 0) when 
the central nucleus is absent, there is no chaos whatsoever, and more than 70% of 
the phase space is dominated by initial conditions of box orbits. However, as mass 
begins to be transferred from the disk to the core, box orbits are depleted and this 
transportation starts to trigger chaotic motion. In particular, we see that as the 
nucleus becomes more and more massive with time the percentage of box orbits 
is significantly reduced and for Mn > 300 it saturates around 25%, while at the 
same time the rate of chaotic orbits exhibits a considerable growth and for about 
t > 6000 time units chaotic orbits is the most populated family, although at high 
enough values of Mn their rate displays a minor decrease. At the end of the galactic 
evolution at t = 10000 time units, or in other words when Mn = 500, the rates of 
box and chaotic orbits tend to a common value (around 25%), thus sharing half 
of the phase plane. Our numerical analysis suggests that the percentages of the 
rest of orbital families change very little during the galactic evolution. Indeed it is 
seen that the 2:1 resonant orbits (meridional bananas) are almost unperturbed by 
the shifting of the mass of the central nucleus occupying about 22% of the phase 
space throughout. Furthermore, the rates of the 1:1 and 4:3 resonant families 
display a monotone behavior roughly around 10% during the mass transportation. 
In addition, we may say that, in general terms, all the other resonant families (i.e.. 


^ It should be pointed out that the color-coded grids of Fig. 4(a-f) show only the A > 0 
part of the phase plane; the A < 0 is symmetrical. Therefore, in many resonances not all the 
corresponding stability islands are present (e.g., for the 2:3 and 4:3 resonances only two (actually 
one and a half) of the three islands are shown). 
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Fig. 5. Time-evolution of the percentages of the different types of orbits in the phase 
(R, R) space of our galaxy model, as a function of the final mass of the nucleus Mn. 


the 2:3, 6:5 and “other”) possess throughout very low percentages (always less than 
5%) thus, the growth of the mass of the nucleus only shuffles the orbital content 
among them. Therefore, taking into account all the above-mentioned analysis, we 
may conclude that in the phase (i?, R) space the types of orbits that are mostly 
influenced by the mass transportation are the box and chaotic orbits. It is also 
interesting to note that the time-evolution of the percentages of the orbits is very 
similar to that given in Fig. 6 in Zotos & Carpintero (2013), where the influence 
of the mass of the nucleus in the time-independent version of same model was 
investigated. 

All the previous numerical results correspond to specific time points off the 
galactic evolution. As it was explained in the previous Section, we can pick some 
characteristic time points, use the respective values of the mass of the disk, the 
mass of the nucleus and the value of the Hamiltonian and integrate orbits in the 
time-independent frozen model. Now we would like to investigate the character of 
orbits and the transition from regularity to chaos and vice versa during the mass 
transportation. For this purpose, we use the full time-dependent Hamiltonian as 
follows: we define in the phase space a grid of 50000 initial conditions of orbits 
with Eq = —480 and we vary the mass of the nucleus form Mn(to = 0) = 0 to 
Arn(tfinai = 10000) = 500 recording at each time step of the numerical integration 
the value of SALT Our computations suggest that the nature of the orbits can 
change either from ordered to chaotic and vice versa or not change at all, as mass 
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Fig. 6. Time-evolution of SALI for six different types of orbits in the time-dependent 
Hamiltonian. Details are given in the text. 


is transported from the disk and a massive nucleus is developed in the central 
region of the galaxy. In Fig. 6(a-f) we present the time-evolution of SALI of six 
different orbits, as the total mass distribution of the galactic system changes with 
time, following the set of equations (3). For all six orbits we set zq = 0, i?o = 0, 
while zq is obtained from the initial value of the Hamiltonian which is Eq = —480. 
The horizontal, blue, dashed line in Fig. 6 corresponds to the threshold value 
(SALI = 10“^) which separates regular from chaotic motion. Fig. 6a shows the 
time-evolution of SALI for an orbit with Rq = 7.5. It is seen that the particular 
orbit starts as regular and remains regular throughout the galactic evolution. The 
pattern of SALI shown in Fig. 6b on the other hand, suggests that the orbit with 
Ro = 0.145 maintains its chaotic character regardless the mass transportation. 
In Fig. 6c where Rq = 14.35, we see that the orbit starts as regular but after 
about 3800 time units it becomes chaotic. The exact opposite scenario is shown in 
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Fig. 7. The shape of the orbit discussed in Fig. 6c in four time intervals of the galactic 
evolution, (a-upper left): 0 < t < 500; (b-upper right): 3500 <t< 3700; (c-lower left): 
3700 < t < 4100; (d-lower right): 8000 < t < 8500. 

Fig. 6d where Ro = 0.23. This orbit exhibits chaotic behavior for the first about 
4800 time units of the galactic evolution however, for t > 4800 it clearly becomes 
regular. The transition from regularity to chaos and vice versa may occur more 
than one time during the mass transportation. Fig. 6e shows such a characteristic 
example of an orbit with Rq = 0.2217 which display s chaotic nature only in the 
interval 1950 < t < 5670. Nevertheless, the determination of the character of an 
orbit is not always very easy. This becomes evident by inspecting Fig. 6f where 
one may observe that the SALI of the orbit with Rq = 0.21 oscillates around 
the threshold value, thus preventing us from having a clear and definitive view 
regarding the character of this orbit, which probably remains sticky throughout 
the galactic evolution. We must point out that these dynamical transitions are 
not related by no means to stickiness or ordinary diffusion phenomena that occur 
in time-independent systems. 

In order to have a more enlightening picture about the transition from order to 
chaos during the mass transportation we provide in Fig. 7(a-d) the shape on the 
meridional (i7, z) plane of the orbit explained in Fig. 6c for four time intervals of 
the galactic evolution. Fig. 7a shows the orbit for the first 500 time units, where it 
is clearly seen that the orbit is beyond any doubt a regular 6:7 resonant orbit. In 
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Fig. 8. Orbital structure of the {R, Mn)-plane. This diagram gives us a detailed 
analysis of the evolution of orbits starting perpendicularly to the i?-axis when the mass 
of the nucleus varies in the interval Mn G [0, 500] during the galactic evolution for 0 < 
t < 10000 time units. 

Fig. 7b where 3500 <t< 3700 we observe that the path of the trajectory the star 
follows becomes very unclear which is an indication of imminent chaotic motion. 
Indeed in Fig. 7c where 3700 < t < 4100 we see that the transition from regularity 
to chaos has been completed, according to the corresponding time-evolution of 
SALI which reported the transition point around 3800 time units. Moreover, it is 
observed in Fig. 7c that the test particle (star) spends a great deal of time moving 
very close to the galactic plane. Finally, in Fig. 7d where 8000 <t< 8500, that is 
an advanced stage of the galactic evolution where the nucleus is massive enough, 
the complete chaotic nature of the orbits is fully revealed and the star moves at 
relatively high distances from the galactic plane up to about 10 kpc. It should 
be stressed out however, that the shape of an orbit gives only fast and qualitative 
information which sometimes can be inconclusive or even misleading regarding the 
nature of the orbit. Therefore, only highly accurate methods that use certain and 
objective numerical criteria, such as the SALI, should be used for determining the 
characters of orbits. 

The color-coded grids of initial conditions in the phase (i?, R) plane presented 
in Fig. 4(a-f) provide information on the phase space mixing for only a specific 
time point of the galactic evolution and for the corresponding value of the mass of 
the nucleus Mn. Following Henon’s idea Henon (1969) however, we can consider a 
plane which provides information about regions of regularity and regions of chaos 
using the section z = R = 0, i > 0, i.e., the test particles (stars) are launched on 
the i?-axis, parallel to the z-axis and in the positive 2 ;-direction. Thus, in contrast 
to the previously discussed grids (Fig. 4(a-f)), only orbits with pericenters on the 
i?-axis are included and, therefore, the value of Mn is now used as an ordinate. In 
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this way, we can monitor how the mass of the nucleus influences the overall orbital 
structure of our dynamical system using a continuous spectrum of values of 
rather than a few discrete ones. Fig. refRt shows the orbital structure of the 
{R, Mn)-plane, when Mn G [0, 500] and 0 < t < 10000 time units. In order to be 
able to monitor with sufficient accuracy and details the evolution of the families of 
orbits, we defined a dense grid of 10^ initial conditions in the (i?, Mn)-plane. For 
creating this plane we use the time-independent model with the values of Mn, Md 
and E according to Eqs. (3) and (9). It is evident, that the vast majority of the grid 
is covered either by box or 2:1 resonant orbits, while initial conditions of chaotic 
orbits are mainly situated to right outer part of the (i?, Mn)-plane. Furthermore, 
the 2:3, 4:3 and 6:5 resonances produce thin vertical stability layers. Furthermore, 
our numerical calculations indicate that in the interval 396 < Mn < 407 there 
is no indication of the 2:3 resonance. This justifies why the stability islands of 
this resonance were found absent in Fig. 4e when Mn = 400. It is also observed, 
that several families of higher resonant orbits are present, corresponding to thin 
filaments of initial conditions living inside the box region. We would like to note, 
that the maximum value of the R coordinate (i?max) is slightly reduced as the 
nucleus gains mass. We must also point out that the (i?, Mn)-plane contains only 
such orbits starting perpendicularly to the i?-axis, while all types of orbits whose 
initial conditions are pairs of position-velocity (i.e., the 1:1 resonant family) are 
obviously not included. 

5. CONCLUDING REMARKS 


In the present work we have sought to shed some light on the interesting 
phenomenon of mass transportation by investigating the orbital dynamics of a 
mean filed galaxy model, when the mass parameters are linearly changing in time. 
Eor this purpose, we used an analytic, axially symmetric, time-dependent galactic 
gravitational model which embraces the general features of a disk galaxy with 
a dense, massive, central nucleus. During the galactic evolution the total mass 
of the galaxy remains constant which means that whatever mass the disk loses, 
it is gained by the nucleus. In order to simplify our numerical calculations we 
chose to work in the meridional (R, z) plane, thus reducing three-dimensional to 
two-dimensional motion. We kept the values of all the other parameters constant, 
because our main objective was to determine the influence of the mass of the 
nucleus on the percentages of the orbits, where mass is transported from the disk 
to the nucleus. Our thorough and detailed numerical analysis suggests that the 
level of chaos, as well as the different regular families, are indeed very dependent 
on the galactic evolution. Eurthermore, we wanted to prove that transitions from 
regularity to chaoticity and vice versa are possible in this simple model. 

Since a distribution function of the model potential was not available so as to 
use it for extracting different samples of orbits, we had to follow an alternative 
path. We defined for several time points of the galactic evolution, dense grids of 
initial conditions (Rq^Rq) regularly distributed in the area allowed by the corre¬ 
sponding value of energy on the phase space. To show how the mass transportation 
influences the orbital structure of the system, we presented for each case the color- 
coded grids of initial conditions, which allow us to visualize what types of orbits 
occupy specific areas in the phase space. Each orbit was numerically integrated 
in the time-independent Hamiltonian for a time period of 10^ time units (10^^ 
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yr), which corresponds to a time span of the order of hundreds of orbital periods. 
The particular choice of the total integration time was made in order to eliminate 
sticky orbits (classifying them correctly as chaotic orbits) with a stickiness at least 
of 100 Hubble times. Then, we made a step further in an attempt to distribute 
all regular orbits into different families. Therefore, once an orbit has been char¬ 
acterized as regular applying the SALI method, we then further classified it using 
a frequency analysis method. For the numerical integration of the grids with the 
initial conditions of the orbits, we needed about between 5 and 7 days of CPU 
time on a Pentium Dual-Core 2.2 GHz PC, depending on the rate of regular orbits 
in each case. 

The most important outcomes of our numerical investigation can be summa¬ 
rized as follows: 

• Numerous types of ordered orbits were identified in our disk galaxy model, 
while there are also extended chaotic domains separating the areas of regular 
motion. In particular, a plethora of resonant orbits (i.e., 1:1, 2:1, 2:3, 4:3, 6:5 and 
other resonant orbits) are present, thereby enriching the orbital structure of the 
galaxy. It should be clarified that by the term “other resonant orbits” we refer 
to resonant orbits with a rational quotient of frequencies made from integers > 5, 
which of course do not belong to the main families. 

• It was observed that the galactic evolution, where mass is linearly transported 
from the disk to the central nucleus, influences mainly the percentages of box and 
chaotic orbits in the phase (i?, R) space. The mass of the central nucleus, although 
spherically symmetrical and therefore maintaining the axial symmetry of the entire 
galaxy, can generate substantial chaotic phenomena in the meridional plane, as it 
is above zero. 

• We found that as the galactic nucleus becomes more and more massive the 
percentage of chaotic motion grows mainly at the expense of box orbits, while 
chaotic orbits is the dominant family once the mass of the nucleus has reached 
about 7% of the mass of the galactic disk. At early stages of the mass transporta¬ 
tion, where the mass of the nucleus is still relatively low, we measured the largest 
amount of ordered orbits. In fact, when <100 more than half of the phase 
space is covered by initial conditions of box orbits. 

• Our results strongly indicate that in the time-dependent Hamiltonian system, 
where a massive nucleus is developed in the central region of the disk galaxy 
through the mass transportation, the character of the orbits can change either 
from ordered to chaotic and vice versa or not change at all. 

• The SALI method was proved highly efficient and accurate in the identifica¬ 
tion of chaos in the time-dependent system. Our computations revealed that this 
indicator is especially suited for detecting time intervals where an orbit exhibits a 
fundamental change in its character. Specifically, by following the time-evolution 
of SALI one can determine in detail the orbit’s successive transitions from regu¬ 
larity to chaoticity and vice versa. 

Judging by the interesting findings we may say that our task has been suc¬ 
cessfully completed. We hope that the present numerical analysis and the cor¬ 
responding results will be useful in the field of time-dependent galactic Hamilto¬ 
nian systems. This is a promising step in the task of understanding the galactic 
evolution of disk galaxies with spherical nuclei. Taking into account that our re¬ 
sults are encouraging, we are planning to properly modify our galactic model in 
order to expand our investigation into three dimensions and explore the entire 
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six-dimensional phase space. In addition, N-body simulations may elaborate the 
mass transportation and its implication on the evolution of galaxies. 

ACKNOWLEDGMENTS. The author would like to express his warmest 
thanks to the anonymous referee for the careful reading of the manuscript and 
for all the apt suggestions and comments which allowed us to improve both the 
quality and the clarity of the paper. 

REEERENCES 

Allen C., Santillan A. 1991, RevMexAA, 22, 255 
Binney J., Spergel D. 1982, ApJ, 252, 308 
Binney J., Spergel D. 1984, MNRAS, 206, 159 

Binney J. Tremaine S. 2008, Galactic Dynamics, Princeton Univ. Press, Princeton, 
USA 

Garanicolas N.D., Papadopoulos N.J. 2003, A&A, 399, 957 
Garanicolas N.D., Zotos E.E. 2013, PAS A, 30, 49 
Garpintero D.D., Aguilar L.A. 1998, MNRAS, 298, 1 

Garpintero D.D., Maffione N., Darriba L. 2014, Astronomy and Gomputing, 5, 19 
Gollin S., Zahn J.P. 1999, A&A 344, 449 

Gomez E., Helmi A., Brown A.G.A., Li Y.S. 2010, MNRAS, 408, 935 

Hasan H., Norman G.A. 1990, ApJ, 361, 69 

Hasan H., Pfenniger D., Norman G.A. 1993, ApJ, 409, 91 

Henon M. 1969, A&A, 1, 223 

Irrgang A., Wilcox B., Tucker E., Schiefelbein L. 2013, A&A, 549, A137 
Manos T., Bountis T., Skokos, Gh. 2013, Journal of Physics A, 46, 254017 
Manos T., Machado, R.E.G. 2014, MNRAS, 438, 2201 
Merritt D. 1999, PASP, 111, 129 
Miyamoto W., Nagai R. 1975, PASJ, 27, 533 
Muzzio J.G., Garpintero D.D., Wachlin E.G. 2005, GeMDA, 91, 173 
Ollongren A. 1962, Bulletin of the Astronomical Institutes of the Netherlands, 16, 
241 

Papadopoulos N.J., Garanicolas N.D. 2006, Baltic Astronomy, 15, 561 
Press H.P., Teukolsky S.A, Vetterling W.T., Elannery B.P. 1992, Numerical 
Recipes in EORTRAN 77, 2nd Ed., Gambridge Univ. Press, Gambridge, USA 
Sidlichovsky M., Nesvorny D. 1996, GeMDA, 65, 137 
Skokos C. 2001, J. Phys. A Math. Gen., 34, 10029 
Zotos E.E. 2011, New Astronomy, 16, 391 
Zotos E.E. 2012, New Astronomy, 17, 576 
Zotos E.E. 2013, Nonlinear Dynamics, 73, 931 
Zotos E.E. 2014a, A&A, 563, A19 
Zotos E.E. 2014b, Baltic Astronomy, 23, 37 
Zotos E.E., Garpintero D.D. 2013, GeMDA, 116, 417 
Zotos E.E., Garanicolas N.D. 2013, A&A, 560, AllO 
Zotos E.E., Garanicolas N.D. 2014, Nonlinear Dynamics, 76, 323 



